###nomo
#引用包
library(survival)
library(regplot)
library(rms)
riskFile="rs_tcga.txt"      #风险文件
cliFile="clinical.txt"       #临床数据文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\8.NOMO\\2")     #设置工作目录

#读取风险输入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)

#读取临床数据文件
cli=read.table(cliFile, header=T, sep="\t", check.names=F)
cli=cli[!duplicated(cli$id),]
row.names(cli)=cli$id
cli=cli[,-1]
cli=cli[apply(cli,1,function(x)any(is.na(match('unknow',x)))),,drop=F]
#cli$Age=as.numeric(cli$Age)

#合并数据
samSample=intersect(row.names(risk), row.names(cli))
risk1=risk[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(risk1[,c("futime", "fustat", "riskScore")], cli[,-c(1,2)])
#rt=rt[,colnames(rt)[c(1:2,4:6,3)]]
rt=rt[order(rt$riskScore),]
rt=na.omit(rt)
#rt=rt[-c(nrow(rt),(nrow(rt)-1)),]
#绘制列线图
#rt=rt[,-5]
res.cox=coxph(Surv(futime, fustat) ~ . , data = rt)
#?regplot
nom1=regplot(res.cox,
             plots = c("spikes", "spikes"),
             spkcol="#78302E",
             clickable=F,
             title="",
             points=TRUE,
             droplines=TRUE,
             observation=rt[2,],
             rank="sd",
             failtime = c(1,3,5),
             prfail = F)

#列线图的风险打分
nomoRisk=predict(res.cox, data=rt, type="risk")
risk1=risk1[rownames(rt),]
rt=cbind(risk1, Nomogram=nomoRisk)
outTab=rbind(ID=colnames(rt), rt)
write.table(outTab, file="nomoRisk.txt", sep="\t", col.names=F, quote=F)

#校准曲线
pdf(file="calibration.pdf", width=5, height=5)
#1年校准曲线
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=1)
cal <- calibrate(f, cmethod="KM", method="boot", u=1, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1),
     xlab="Nomogram-predicted OS (%)", ylab="Observed OS (%)", lwd=3, col="#E64B35FF", sub=F)
#3年校准曲线
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=3)
cal <- calibrate(f, cmethod="KM", method="boot", u=2, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", lwd=3, col="#4DBBD5FF", sub=F, add=T)
#5年校准曲线
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=5)
cal <- calibrate(f, cmethod="KM", method="boot", u=3, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="",  lwd=3, col="#00A087FF", sub=F, add=T)
legend('bottomright', c('1-year', '3-year', '5-year'),
       col=ggsci::pal_npg("nrc")(10)[1:3], lwd=3, bty = 'n')
dev.off()

